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Grease ice is an agglomeration of disc-shaped ice crystals, named frazil ice, which forms 
in turbulent waters of the Polar Oceans and in rivers as well. It has been recognized 
that the properties of grease ice to damp surface gravity waves could be explained in 
terms of the effective viscosity of the ice slurry. This paper is devoted to the study of 
the dynamics of a suspension of disc-shaped particles in a gravity wave held. For dilute 
suspensions, depending on the strength and frequency of the external wave flow, two 
orientation regimes of the particles are predicted: a preferential orientation regime with 
the particles rotating in coherent fashion with the wave held, and a random orientation 
regime in which the particles oscillate around their initial orientation while diffusing 
under the effect of Brownian motion. For both motion regimes, the effective viscosity has 
been derived as a function of the wave frequency, wave amplitude and aspect ratio of the 
particles. Model predictions have been compared with wave attenuation data in frazil ice 
layers grown in wave tanks. 



1. Introduction 

Grease ice is a thin slurry of disc-like platelets of ice crystals, called frazil ice, which 
forms in supercooled waters of the Polar Oceans under cold and windy conditions. Frazil 
discs measure approximately 0.1 — 0.4cm in diameter and 1 — 100/im in thickness (Kivisild 
1970). Grease ice can accumulate up to tens of centimeters and significantly affect ocean 
surface roughness by attenuating short waves. This effect has been widely documented by 
observations of early whalers. Furthermore, synthetic aperture radar imagery of grease 
ice scenes appears dark because of the suppression of the gravity-capillary waves resonant 
with the incident microwave radiation (1 — 10cm) (Wadhams & Holt 1991). 

Laboratory measurements of wave propagation in grease ice show that wave attenu- 
ation can be explained in terms of the medium effective viscosity (Newyear & Martin 
1997). Martin & Kauffman (1981) developed a viscous-plastic model to explain the ob- 
served wave attenuation. They claimed that the viscous nature of grease ice could arise 
from interactions among frazil crystals leading to the presence of an energy sink in the 
wave dynamics. The authors did not present any estimate of grease ice effective viscosity 
from their data. On the other hand, wave dispersion and attenuation data of Newyear 
& Martin (1997) were consistent with a constant viscosity value, comparable to that of 
glycerin at 0°C, in the range of frequencies from 6.6 to 9.5s -1 (Newyear & Martin 1999). 
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The viscosity was estimated using a two-layer wave propagation model, which represents 
grease ice as a viscous fluid superimposed on inviscid water (Keller 1998). 

As a matter of fact, the concept of bulk viscosity for grease ice holds because the size 
of the frazil particles (~ 0.1cm) is much smaller than both the vertical scale of grease 
ice (10cm in the laboratory) and the horizontal scale (from ~ 100cm in the laboratory 
to hundreds of meters in the ocean) of the travelling wave (Newyear & Martin 1999). 

From the theoretical point of view, it is possible to have detailed information on the 
behavior of a suspension of disk-like particles in the velocity field of a gravity wave only in 
the dilute limit. In this limit, the problem becomes that of the behavior of an individual 
particle in a given flow field, in the absence of interactions with the other particles in 
suspension. Iterative approaches such as that of the differential scheme (Bruggeman 1935) 
can then be used to obtain semi-quantitative informations in the high volume fraction 
regimes characteristic of grease ice. 

A second simplification is obtained disregarding inertia effects at the scale of the ice 
platelets. Starting from the work of Jeffery (1922), much work has been devoted to the 
dynamics of an ellipsoidal particle under creeping flow conditions. The importance of 
Brownian motion for the presence of an equilibrium particle orientation distribution, 
and consequently for the existence of a uniquely defined bulk viscosity was already rec- 
ognized in (Taylor 1923). It turns out that, unless the particles are so small that Brownian 
diffusion dominates, the external flow strongly influences the orientation distribution and 
this effect cannot be disregarded for relatively large particles like the ice platelets. As rec- 
ognized in (Bretherton 1962), this may lead, in flow regions characterized by high strain 
and low vorticity, to the possibility of fixed orientation regimes for the particles. Only 
more recently, however, have time dependent situations involving ellipsoids in suspen- 
sion, come under scrutiny. In (Zhang & Stone 1998), the forces and torques acting on an 
oscillating disk in a quiescent fluid have been calculated. In (Szeri et al. 1992), the ori- 
entation dynamics of an ellipsoidal particle under the effect of combined time-dependent 
vorticity and strain has been analyzed. 

The calculation of the bulk viscosity of a dilute suspension of ellipsoids in a plane 
shear was carried on in (Leal & Hinch 1972), in the small but non-zero Brownian motion 
regime. The effective viscosity of a concentrated suspension of aligned disks was studied in 
(Sundararajakumar et al. 1994), using slender body theory arguments. In (Phan-Thien 
& Pham 2000), a differential scheme approach was used to calculate the effective viscosity 
in the concentrated regime assuming random orientation of the ellipsoids. In all cases a 
time independent situation was considered. 

In the present paper, we shall consider the case of gravity waves in infinitely deep 
water. In this case, it is possible to pass to a reference frame in which the flow is time 
independent, and this eliminates the possibility of irregular orbits in orientation space, 
observed in (Szeri et al. 1992) already in the case of simple periodic flows. 

This paper is organized as follows. In the next section, the orientation dynamics of 
a Stokesian particle in a deep water gravity wave will be elucidated. In particular, the 
possibility of coherent collective motions in the suspension will be examined. In Section 
III, the effective viscosity of a dilute ellipsoid suspension will be calculated, analyzing its 
dependence on the wave frequency and amplitude. In section IV the merit and limitation 
of a creeping flow approach to modelling the frazil ice dynamics will be discussed. In 
Section V, using a differential scheme, the results will be extrapolated away from the 
dilute limit, and will be compared with available data from wave-tank experiments. 
Section VI will be devoted to conclusions. 
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2. Orientation of a disk-like particle in the velocity field of a deep 
water gravity wave 

We consider a dilute suspension of rigid oblate axisymmetric ellipsoids, supposed small 
enough that inertia be negligible at the particle scale. We also suppose that the particles 
are free of external forces or torques and that the suspension is so dilute that the effects 
of mutual interaction among particles are negligible. In this dilute limit, the rheological 
properties of the suspension will descend from the response of a single particle to the 
time dependent wave flow. 

In order to represent the wave flow, we introduce a reference frame with the origin 
at the water surface, the Xi-axis along the direction of propagation of the wave and the 
a^-axis pointing vertically towards the sea-bottom. In the hypothesis of small amplitude 
inviscid waves in infinitely deep water, we obtain the following velocity field: 

U\ = U exp( — kx2) sin(fca;i — ojt) ,^ ^ 

U-2 = U exp(— kx2) cos(fcxi — ujt) 

where U is a typical value for the fluid velocity in the wave field. 

The orientation dynamics in the presence of fore-aft symmetry, will in turn be de- 
termined by the balance between the strain rate E and the vorticity ft of the wave at 
the particle position, again provided the particle is sufficiently small to allow lineariza- 
tion of the wave field on its scale. In the case of a revolution ellipsoid, with symmetry 
axis identified by the versor p, the orientation dynamics will be described by the Jeffery 
equation: 

p = n p + G[E • p (p • E p)p] + 0{{kaf) (2.2) 

where G is the ellipsoid eccentricity defined in terms of the particle aspect ratio r = a/6, 
where a and 6 are respectively along and perpendicular to the symmetry axis, by means 
of the relation 

G = -j— ^ 

r 2 + 1 

For disk-like particle we have clearly r <C 1 and G ~ — 1. For small amplitude waves, 
the particle displacement will be small with respect to A: -1 and we can approximate the 
instantaneous value of the strain felt by the particle, with its value measured at the 
initial position x. Furthermore, for linearized waves, the induced velocity field is confined 
within a region whose thickness is of the order of the wave amplitude (defined as the 
valley to crest semiheight): A <C fc _1 . We can then approximate exp(— kx^) = 1. From 
Eq. H2.1f> we find easily the expression for the strain rate: 

E = ku( cos (^i-^) -sfo(fczi-wi) A {2 3) 

\ — sm(rai — ujt) — cos(kxi — ujt) J 

while the vorticity 1~2 is identically zero thanks to the potential flow nature of the inviscid 
wave field. Equation i|2.3[) describes a strain field rotating with frequency cu/2 around 
the X3 axis. Changing variables to a reference frame rotating with the strain, the time 
dependence in E disappears and a non-zero vorticity is produced: 



(we identify components in the rotating frame with an overbar). For each value of x±, we 
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Figure 1. Orientation of an ellipsoidal particle in a strain field rotating with angular velocity 
to with respect to the laboratory [x\ axis). For ^ u> ^ 1 the symmetry axis p is confined to 
the strain plane x 1 X2- For Co — > 0, alignment with the strain instantaneous compressive direction 
ip = — 7r/4 occurs. For < Co ^ 1, the symmetry axis of the particles lags behind by a constant 
angle V, with ij) = —n/2 for w = 1. For a) > 1 no stationary solution for ?/> exists, corresponding 
to the particle being unable to follow the rotating strain. 



choose the new variables in such a way that the strain rate reads: 

E = kU ( J J ) (2.5) 

corresponding to the strain expansive direction placed at 7r/4 with respect to the rotating 
X\ axis. Introducing polar coordinates (see Fig. [Q, and normalizing time and vorticity 
with the strain strength e = kU, Jeffery's Eq. 1)2. 2f) leads to the following system of 
equations: 

ijj = -£> - cos2V> ,„ g s 

= -|sin26»sin2^ 1 ' ' 

where Co = —io/2Ge and / = d//df, f = —Get. For tD < 1, this system of equations has 
equilibrium solutions (^, 0) = ^(cos -1 Co — nir, rmr). Of these, only the one 

ip = — cos -1 Co — — + nir, 8 = — + nir (2-7) 

is stable and is approached in a time ~ e _1 . For Co > 1, instead, choosing the time so 
that ip(0) — 0, we have the trajectories: 

tamfr(t) = -(§±1) 1 tan[(^ 2 - 1)**], 

tan ^X(- T ^ I y)"tane(0) 

We thus have a high strain regime in which, as illustrated in Fig. ^ the particles in 
suspension are all aligned with the local strain and rotate in coherent fashion, and a low 
strain regime in which the particles do not rotate, rather, they oscillate around their 
natural orientation. As it is easy to see from Eq. JUS}, the transition from the low to 
the high strain regime is characterized by the particle spending an increasing amount of 
time, as Co — > 1, near (— 7r/2, 7r/2). This corresponds to the rotation period of the particle 
(Co 2 — 1)~2 (always measured in the rotating frame) going to infinity, as tp = — n/2 
becomes a fixed point for the system. 

The linearized gravity waves theory allows us to write e = kU in terms of the wave 
amplitude A defined as the valley to crest semiheight, the gravitational acceleration g 
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and the wave frequency u), starting from the dispersion relation: 

UJ 2 

k = — (2.9) 
9 

and the expression of the typical wave velocity value: 

U = Alj (2.10) 

This leads to the expression for the strain strength e in the small r regime corresponding 
to G = — 1 : e = and to the condition for the existence of equilibrium 




(2.11) 



We thus see that in the case of gravity waves, the aligned particle case corresponds to a 
high frequency (or small wavelength) limit. The same transition is observed experimen- 
tally in the case of grease ice (see pages 307 and 308 in Martin & Kauffman 1981), with 
the crossover frequency at uj ~ ^/0.35.g/A 



3. The bulk stress and the effective viscosity of the fluid-particle 
mixture 

The bulk stress of a dilute suspension of axisymmetric ellipsoidal particles is given by 
the law (Leal & Hinch 1972), (Hinch & Leal 1975): 

<r = PI + 2^E + 2^{2A(pppp) : E + 2B[<pp) • E + E • (pp)] +CE + F(pp) • D} (3.1) 

where /i is the dynamic viscosity of the pure fluid, P is the pressure, (f> is the volume 
fraction of the particles, A,B,C,F are dimensionless shape coefficients and D is a term 
that takes into account Brownian motion effects. It is an open question whether other 
effects, such as interaction with other particles, could be modelled by a noise term. The 
presence of this term, independently of its amplitude, guarantees that memory of any 
initial particle orientation, including unstable equilibrium points, is lost and a statistical 
equilibrium state, in an 0(D~ r ) time, is eventually reached. 

Following (Leal & Hinch 1972), we shall consider the small noise limit in which D^ 1 is 
much longer than the other timescales of the process, which are given in dimensionless 
form by (a) 2 — . Over these timescales, the evolution of the process will be therefore, 
to lowest order, that of the unperturbed system. 

The second and fourth moment of p are calculated function of the PDF (probability 
density function) for the particle orientation p(0,ip,t). The A,B,C coefficients may be 
obtained from (Jeffery 1922), in terms of the following elliptic integrals: 

dX f°° XdX 



a = 



and 



o (6 2 +A) 3 V7+A' Jo (b 2 + XfVoJTX 

XdX 



ff = r - F=f 

7o (P + X)^a 2 + X)Va T TX , J 



(b 2 + X) 2 (a 2 + X)Va T Tx' Jo (b 2 + X) 2 (a 2 + X)VoJTX 

where a and b identify the ellipsoid semiaxes parallel and perpendicular, respectively, to 
the symmetry axis. More precisely 

a" 1 2 _ 1 1 , _ 1 

"F TTTo 7 7T77 o ', TtT: B = -——-r- - —— and C 



2b 2 a>f3" 2b 2 a> /?'(a 2 F6 2 )' /3'(a 2 + b 2 ) b 2 a' b 2 a' 
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In the case of disk-like (r <C 1) particles, disregarding 0(r) terms: 

5 104 „ 4 64 1 , „ 8 128 

A= 3^ + 9^- 1 ' B = - 3^-9^ + 2 and C =3^ + 9^ < 3 ' 2 > 

Notice that this value of C differs at subleading order 0(1) from the one in (Leal & Hindi 
1972). As seen in the previous section, two orientation dynamics regimes are possible and 
these affect the value of the angular averages entering Eq. $A.1\ . We consider in detail 
the two regimes below. 

From the stress <x, it is possible to calculate an effective viscosity ft, in terms of the 
viscous dissipation in the suspension, exactly as it is done with spherical particles: 

fi=l^^~0- + K4)iM (3.3) 
where K is called the reduced viscosity for the suspension. 

3.1. Preferential orientation regimes: ^ Cj ^ 1 

In this regime, after a relaxation time ~ e _1 , all particles tend to align, in the rotating 
frame, in the direction identified by Eq. 1)2. 7|) . For small diffusivities, the variance of 
the distribution around these fixed points will be D/e. As already mentioned, this state 
of affairs corresponds, in the laboratory frame, to the particles rotating in coherent 
fashion with the wave field. The fourth and second order tensors (pppp) and (pp) have 
a simpler form in the rotating reference frame with the x\ axis along p. In this new frame 
of reference, the rate of strain tensor E takes the following form 



E 



■VT 



lu yl — lu 2 



while the (pppp) and (pp) tensors read: 

{piPjPkPi} = SiiSijSikSu and (pipj) = hihj 

where Sij is the Kronecker delta. Substituting into Eqs. I|3.1|l and (|3.3|) . the reduced 
viscosity coefficient K, is promptly obtained: 

K = A(l~ui 2 ) + 2B + C (3.4) 

From Eq. (|3.4(l . the dominant 0(r _1 ) contribution to the viscosity is the u> dependent 
contribution proportional to A, while 2B + C = 1 + 0(r). For this reason, the reduced 
viscosity K is characterized by a minimum at the crossover w = 1, at which K ~ 1 
[compare with the spherical particle value K = 5/2, (Landau 1959)]. 

3.2. Continuously rotating regime: Co 1 

In the laboratory frame this regime corresponds to the particle oscillating around its ini- 
tial orientation, while slowly diffusing with respect to angle, under the effect of Brownian 
couples. In the rotating frame, the problem can be mapped to that of the ellipsoid in a 
plane shear: the equation of motion for a particle in the rotating frame (|2.8J) is in fact 
identical to that of a particle with aspect ratio 



(^tt)* (3 ' 5) 

\0J + 1/ 



in a plane shear lu = 2e. The equilibrium distribution of an ensemble of particles whose 
orientation dynamics is described by Eq. I|2.6|l . in the presence of an isotropic Brownian 
couple, is then obtained from the theory of (Leal & Hinch 1972), whose main results are 
reported below. 
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The particle orientation is identified by the variables t and c where t is defined by the 
first of Eq. IZBJl and gives the normalized time needed, on the Jeffery orbit starting from 
the current values of 9 and i/j, to go from -0 = to the current value of tp, while c obeys: 



/(j + cos2?/a 
V oj-1 ) 



tan 9 



Thus, tan 1 c is the value of 9 at tp — 7r/2, and identifies the Jeffery orbit. 
In these variables, the orientation PDF can be decomposed as: 

p(c,i) = p{i\c)p(c) 



(3.6) 



(3.7) 



where, from the fact that t is a time along a trajectory, p(t\c)dtdc is also the infinitesimal 
solid angle element in the variables t and c. The marginal PDF p(c) is given by: 



where 



and 



p(c) = const. c[(iJ 4 c 4 + H 2 c 2 + H )F}- 



s 2 + 1, 



Ho, 



1 

4^2 



Ho 



r(* 2 + l) 



(3.8) 



2H iC 2 +H 2 -S 



2H iC 2 +H 2 +S 
2(g 2 -4) 



H 2 > AH 4 H , 



eX P \_2HiC 2 +H 2 

expos'" 1 (H 2 



where S= \H 2 - AH 4 H () \ t . 



Hi 



AH A H 



(3.9) 



4) tan- 1 S , - 1 (2 J ff 4 c 2 + fl" 2 )] 



H 2 < 4H 4 H , 



The PDF p(c) is actually the only thing that we need, since the averages along the 
orbits of the tensors pp and pppp are already available (Jeffery 1922). In fact, from Eq. 
H3.3(l . the reduced viscosity can be written as 



K = A(sin 4 9 sin 2 20) + 25 (sin 2 6) + C 



(3.10) 



but, from (Jeffery 1922): 



(sin 4 9 sin 2 2tp\c) 



2s 2 



(s 2 - If 



c 2 (.s 2 + l) + 2 
[(c 2 s 2 + l)(c 2 + l)]i 



and 



•9\c) 



[(C 2 S 2 + 1)(C 2 + 1)]2 

Completing the averages by means of Eq. (|3.8|l , leads to behaviors for (sin 4 9 sin 2 2ip) 
and (sin 2 9) shown in Fig. [3 Substituting into Eq. I|3.1U|) with the expressions for the 
coefficients A, B and C provided by Eq. (|3.2|l . allows to determine the reduced viscosity 
of a dilute suspension of ellipsoidal particles, for arbitrary values of the aspect ratio r 
and of the reduced frequency a). As in the continuously rotating regime, we see that the 
effective viscosity grows away from the crossover at 0) = 1, with the dip becoming more 
pronounced as the aspect ratio r is sent to zero. This is illustrated in Fig. [21 in the case 
of a disk-like particle with a value of the aspect ratio in the range characteristic for frazil 
ice. Notice that the asymptotic regime of random particle orientation p(9, ip) 
leading to the expression for the reduced viscosity (Phan-Thien & Pham 2000) 



sin V 

47T 



4 4 
K= — A+-B 
15 3 



C 
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Figure 2. Plot of (sin 4 9 sin 2 2i/>) (bottom) and (sin 2 9) (top) vs. Co in the low strain range. 



K 




o 1 ' 1 1 1 

0.5 1 1.5 2 



Figure 3. Reduced viscosity vs. normalized frequency Co ~ oj/2e for a disk-like particle with 
aspect ratio r — 0.045. From Co oc a; -2 , the long wave regime corresponding to random particle 
orientation, occurs for Co > 1 and the short wave one, corresponding to coherent motion, for 
Co < 1. The horizontal line to the right gives the asymptotic value K(Co — » oo). 

is obtained already for relatively small values of the reduced frequency Co ~ 2. 

We remark that, using the expression for C given in (Leal & Hinch 1972), would have 
produced an unphysical negative value of the reduced viscosity K at the crossover. 

4. Grease ice 

The analysis carried on so far assumed a dilute regime, which is very different from 
the conditions typical of grease ice. Furthermore, the analysis assumed creeping flow 
conditions, which may be problematic for millimiter size particles. As regards the first 
issue, the differential scheme, originally proposed by (Bruggeman 1935), provides an 
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analytical method to generalize the well-known Einstein formula for the effective viscosity 
of dilute suspensions to finite concentrations (Brinkman 1952), (Roscoe 1952). More 
recently, the differential scheme has been exploited to study the viscosity problems related 
to randomly oriented spheroidal inclusions in viscous fluids (Phan-Thien & Pham 2000). 

The basic idea is to calculate the increment of effective viscosity that is obtained adding 
a volume Av of particles to a volume vo of suspension already containing a fraction v/vq 
of particles. Let us indicate with v = p/p w the kinematic viscosity of the solvent and 
with v ~ p,/p w the same quantity referred to the suspension, with p w the density of 
the solvent, assumed approximately equal to that of the suspension. The increment in 
volume fraction is 

A = V + A a ~ — - 0- - v/v )Av/v = (1 - <j))Av/v 
Vq + Av v 

while the increment in effective viscosity v will be 

KvAv 



Av 



v + Av 

leading to the differential equation 

dv/d(f> = Kv/{1 - <t>), 9(0) = v (4.1) 

The differential scheme assumes implicitly that viscosity renormalization is the only 
effect of the particle inclusion, which is strictly true only when the particle orientation 
distribution and consequently the viscosity tensor are isotropic. In the opposite limit 
to ^ 1, the theory of (Sundararajakumar et al. 1994) could be applied. 

In the case of microscopic (Stokesian) particles, Eq. (|4.1|) could be integrated from the 
initial condition (corresponding to pure solvent) to the final concentration cj>, using for 
K the constant value obtained from Eqs. I|3.2I3.4I3.10[I and the data in Fig. [21 I n this 
case, the solution would be 

P(4>)=u(l-^)- K , (4.2) 

(This equation allows a maximum packing fraction 4>max = 1 which is above the real 
value 0MAX — tt/4 appropriate for stacked discs). In the case of frazil particles, in 
the integration of Eq. (|4.1|l . there is an initial range of values of <j) in which P(4>) is 
likely to be too small for the creeping flow approximation to apply; in that range, the 
stress coefficients and K will likely differ from their zero Reynolds number limit. Hence 
K = K((f) and Eq. (|4.2|l will provide at most an order of magnitude estimate for the 
effective viscosity of grease ice. 

To determine the importance of these effects, let us consider an individual ice platelet 
in a suspension of effective viscosity D. Creeping flow conditions require stationarity and 
small particle Reynolds numbers. Introducing the effective Stokes time fs ~ b 2 /v: 

ujb 2 eb 2 
cot s = < 1, Re p = — = ujT S kA < 1, (4.3) 

V V 

where e = kU is the strain strength. In the case of a dilute suspension, v = v ~ 
0.01cm 2 s~ 1 ; taking for the particle radius the value b ~ 0.1cm, we would obtain from 
Eqs. H2.9I2.1U[) . values for lots ranging from 1 in open sea, to 10 in laboratory conditions. 
As regards the condition on Re p , this is satisfied in open sea, where kA is of the order 
of a few hundredths, but only marginally in the laboratory, where kA ~ 0.3 (Martin & 
Kauffman 1981). Clearly, the conditions of Eq. (|4.3|l are going to be satisfied in the case 
the suspending medium is the grease ice, when p, ~ 10 2 cm 2 s _1 . This has the consequence, 
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in particular, that the results on the orientation dynamics of Section II are expected to 
remain overall valid. 

In the dilute case, the non satisfaction of the conditions in Eq. (|4.H|> has another conse- 
quence, namely, inertia will cause relative particle-fluid motions and additional dissipa- 
tion in the suspension. We can obtain an estimate of this effect. The particle velocity V, 
subtracted of the contribution from the buoyancy produced drift, will obey an equation 
in the form: 

d 1 \, , dV 1 , , 

— + — n- (V-U)-er— + — F-U + f 4.4 
\at ts ) at ts 

with e = 1 — Pp/pw — 0.1, p p indicating the ice density, F = 0((kb) 2 ) accounting for 

the Faxen force and II the adimensionalized resistance tensor, whose components are 

O(l) for the range of ufs we are interested in (Zhang & Stone 1998). The term f is 

a noise contribution accounting for collision effects with other particles and Brownian 

motion. All inertia effects in the wave flow and from the particle relative motion are in 

first instance neglected. 

The contribution to relative motion from collisions is important in bubbly flows (Kang 

ct Al. 1997), but can be shown to be negligible in the present case due to the regime 

Re p < 1. Let us show this. We can estimate the noise amplitude in a kinetic approach 

by introducing first the collision frequency t^ 1 

r" 1 ~ nAVb 2 = 

a 

where n = 4>/ab 2 is the numerical density of the particles, AV ~ |V — U| estimates the 
typical collision velocity and b 2 estimates the collision cross section. Collisions will be 
important provided r c < fs; in this case, taking the noise as uncorrelated (f(t)f(t')) ~ 
DS(t — t'), we would have for its amplitude: 

From Eq. Q4.4|) . considering / dominant on the other terms to right hand side, we ob- 
tain the estimate for the velocity AV ~ e(Dfs)' 1 ; this is the velocity difference in the 
wave field sampled by the diffusing particle in the relaxation time fs- Substituting the 
expressions for D and fs, we find 

AV~°^ 



and substituting in the expression for r c and comparing with fs, we see that the condition 
fg > t c is equivalent to Re p > 1. Interparticle collisions can then be neglected. 

Passing to the contribution to V — U from direct acceleration by the wave field, we 
see that, for waves in both laboratory (k ~ 10 _1 cm _1 ) and open sea (fc ~ 10 _3 cm _1 ) 
conditions, the Faxen force can be disregarded. Neglecting the noise in Eq. (|4.4|) . we 
obtain in this limit: 

|V — U| ~ erll min(l, cofs) 
The dissipation produced by a single particle due to its translation relative to the fluid 
can be estimated from the product of the drag force in Eq. I|4.4|l and | V— U| as p w b 3 H\\~ — 
U| 2 jfs- The contribution from relative particle fluid rotation will be smaller by a factor 
kb. The dissipation per unit volume will be therefore: 

■eU 



W tr ~ p4>ry—j min(l,wf s )' 
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to be compared with the viscous dissipation W v — jl(kU) 2 . We thus see that the contri- 
bution to dissipation from relative particle-fluid motion is dominant in the dilute case, 
but decreases with 1/p, when u>fg < 1. Taking for fi values in the range of the hundreds, 
and 4> ~ 0.3, as observed in grease ice (see also Table 1), we see that the contribution 
to dissipation from relative particle-fluid motion is negligible, with the ratio Wtr/W^ 
varying from 1(T 4 in the open sea to 10 8 in the laboratory. 

In conclusion, inertia and non-stationarity are likely to contribute to the effective 
value of the grease ice effective viscosity, but do not affect the frazil orientation dynamics 
described in Section II. It is also confirmed that the dominant contribution to dissipation 
and to the effective viscosity is the particle induced stress and not particle-fluid motions. 



5. Comparison with experiments 

We have compared the order of magnitude estimate for the effective viscosity, provided 
by Eq. I|4.2I) . with the wave tank data of (Martin & Kauffman 1981). In their experiment, 
concentrated suspensions of grease ice with thicknesses varying from 7 to 15cm and 
volume fraction <j) between .28 to .44, were allowed to grow in a 2m long tank previously 
filled with saline water to a depth of 41cm. We selected those measurements relevant to 
propagation of deep waves in grease ice layer according to the criterion kh ^ 7r/2, where 
k is the open water wavenumber and h is the ice layer thickness (Phillips 1966). The 
relevant parameters of the experimental data we are considering are listed in Table 1. As 
already discussed, the differential scheme assumes an isotropic suspension. Comparing 
the values of to in Table 1 with Figure |3 we see that the data fall in the range where this 
assumption holds. 



to 




<, 


i 


^(s- 1 ) 


.A(cm) 


h{cm) 


(Z(m-i) 


v(cm x s 2 ) 


K 


1.4-=- 


1.5 


0.30- 


-0.34 


14.9 


1.45-=- 


1.55 


6-^ 


7 


5.3 ±0.4 


243- 


-250 


18.7 ±0.2 


1.3-=- 


1.4 


0.29- 


-0.34 


14.9 


1.6-=- 


1.7 


7-=- 


8 


6.6 ±0.6 


353- 


-361 


19.7±0.3 


1.2-=- 


1.3 


0.28- 


-0.32 


14.9 


1.7-=- 


1.8 


8-=- 


9 


5.8 ±0.5 


275 - 


-280 


20.3 ±0.2 


1.2-=- 


1.5 


0.28- 


-0.34 


14.9 


1.5-=- 


1.8 


7 + 1 


1.5 


1.6±0.1 


587- 


-603 


15.6 ±0.1 


1.2-=- 


1.5 


0.34- 


-0.35 


14.9 


1.5-=- 


1.8 


8.5- 


10 


7.5 ±0.5 


454- 


-463 


18.6 ±0.3 


1.4-=- 


1.4 


0.35- 


-0.44 


10.7 


3.0^ 


3.1 


14 -=- 


16 


3.5 ±0.2 


1010- 


- 1029 


17.2 ±0.2 



Table 1 



Grease ice effective viscosities were estimated from the measured spatial attenuation rate 
q = —A~ 1 dA/dxi using a two-layer viscous fluid wave propagation model (De Carolis & 
Dcsidcrio 2002). The model can be inverted to obtain the effective viscosity v of the grease 
ice from the experimentally observed values of the wave frequency u, the attenuation rate 
q and the thickness h and volume fraction <fi of the frazil. The reduced viscosity K can 
then be determined from Eq. (|4.2|l in order to estimate the corresponding particle aspect 
ratio by means of Eqs. I|3.2I3.4I3.10|I and the data in Fig. [21 From the data in Table 1, 
we obtain the estimate: r~2x 10~ 2 , to be compared with the individual observation 
presented in Martin & Kauffman (1981) of a disk diameter of 0.1cm and thickness of 
1 — lO/im (see their Fig. 11). The corresponding range of variability for frazil ice in 
geophysical environment is 0.1 — 0.4cm in diameter and 1-100/im in thickness (Kivisild 
1970), corresponding to 0.25 xl0~ 4 < r < 0.1. 
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6. Conclusions 

We have obtained predictions on the effective viscosity dependence of a suspension of 
disk-like particles, in the velocity field of deep water waves, on the aspect ratio and the 
concentration of the particles. A key parameter appears to be, in the dilute limit, the 
relative strength of the wave field strain, which parameterizes the wave amplitude, and 
the wave frequency. For high amplitude waves, collective alignment of the particles in 
suspension with the wave field is possible, with the crossover to this regime signalled by a 
deep minimum in the effective viscosity. This minimum is actually lower than the effective 
viscosity in the case of spherical particles and can be smaller, by orders of magnitude, 
than the value of the effective viscosity of a disk-like particle suspension away from the 
crossover. 

An interesting question is whether these behaviors are preserved away from the deep 
water wave regime we have considered in this paper. For shallow water waves, a rotating 
system in which the flow becomes time-independent does not exist anymore, and irregular 
behaviors of the kind described in (Szeri et al. 1992) become possible. 

Some of these results extend away from the dilute limit to the case of grease ice, in 
particular, the presence of a crossover to coherent alignment of the particles for large 
amplitude waves (Martin & Kauffman 1981). As regards the effective viscosity of the 
grease ice, this appears to be dominated by the stress contribution from the individual 
particles rather than from momentum transport from relative particle-fluid motion. This 
is in contrast with other situations, e.g. bubble laden flows (Kang et Al. 1997), in which 
the second is the dominant effect. A creeping flow based calculation of the effective 
viscosity, aided by the use of a differential scheme, to deal with the high concentration 
regime, leads to results consistent with experiments. 
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